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Abstract 

We assess the quantitative accuracy of the particle-in-cell (PIC) approximation used 
in recent ab initio predictions of the thermodynamic properties of hexagonal-close- 
packed iron at the conditions of the Earth's inner core. The assessment is made 
by comparing PIC predictions for a range of thermodynamic properties with the 
results of more exact calculations that avoid the PIC approximation. It is shown that 
PIC gives very accurate results for some properties, but that it gives an unreliable 
treatment of anharmonic lattice vibrations. In addition, our assessment does not 
support recent PIC-based predictions that the hexagonal c/a ratio increases strongly 
with increasing temperature, and we point out that this casts doubt on a proposed 
re-interpretation of the elastic anisotropy of the inner core. 
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1 Introduction 



Ab ini t io calculations ba s ed on density functional theory (DFT ) (|Hohenberg and Kohn 
( 1964 ); Kohn and Sham ( 1965| ): Jones and Gunnarssonl (|l989i l) are widely used 



to calculate the properties of materials at the extreme pressures found in 
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the interior of the Earth and other planets (ptixrude et al.l (|1998f )). and are 
known to be capable of high accuracy. For many years, such calculations 
did not generally include thermal effects, and explicitly treated only zero- 
temperature materials. However, the past few years have seen an increasing 
effort to apply DFT calculations to high-temperature solids and liquids of ge- 
ological interest, us i ng ah initio l attice dynamics o r mole c ular-dynam i cs sim - 
ulation JXFe et all feOOll 120021): iBelonoshko et all fcOOOh: ILIio et al.l (1200(1) : 
Stixrude et al.l (r997t):IWasserman et all (|l996l ): lBrodholt et al.l (|2002l ): IOganov et al 



( 20011 ); ISteinle-Neumann et all (J2002J) ) . Because dynamical DFT calculations 
demand large computer resour ces, a simplifying appro x imation known as the 
'particle-in-cell (PIC) method dHirshfelder et al I dl954h:lHolt and Rossi (Il97(t: 
Holt et all (ll970h : lRee and Holtl (|l973h : IWestra and Cowlevl (|l97^ : ICowlev etal 



tool)) has been used in some work on high-T, high-p solids, an important 
example being the very recent use of the PIC appro ximation to re-interpret 



the el a stic a nisotropy of the Earth's solid inner core (iSteinle-Neumann et al 



,HL - — - - ■ - - - - 

challenged (|Alfe et all ( 20011 )) and conclusions based on it are not necessarily 



(|200lL 120021) ). However, the reliability of the PIC method has not gone un- 



secure. To shed light on this question, we present here the results of our 
own PIC calculations of a range of thermodynamic properties of solid iron at 
Earth's core conditions, which we compare with more exact calculations that 
avoid the PIC approximation. 



The purpose of the PIC approximation is to provide a way of calculating 
the free energy of a vibrating crystal. The essence of the approximation is 
that correlations between the vibrational displacements of different atoms are 
neglected, so that each atom is treated as vibrating independently of every 
other, as in the Einstein model for a vibrating solid. Although the approx- 
imation seems at first sight rather crude, empirically it has been shown to 
yield satisfac tory prediction s for t h e thermodynamic properties of a num - 
ber of solids dHolt and Rossi (|l970f ): iHolt et all (|l97flh : iRee and Holtl (|l973l) : 
Westra and Cowlevl $197$i ). at least for temperatur es above the Debye tem- 
perature. In the geological context, it was shown by Wasserman et all (|l996l ) 
that, when implemented using DFT methods, it gives predictions for the p(V) 
relation on the shock Hugoniot in excellent agreement with experiment for 
Fe up to the melting curve. The PIC approach based on DFT has also been 
applied successfully to the h i gh-p, high-T propertie s of ot her transition metals 
(jCulseren and Cohenl (|200ll ): ICohen and Giilsere3 (|2002l) ). 



In the work on the elastic anisotropy of the Earth's inner core (ISteinle-Neumann et al 
(|200lLl2002ln mentioned above, PIC was used to calculate the elastic constants 
of hexagonal-close-packed (hep) Fe over a range of pressures and temperatures. 
A key prediction was that the c/a ratio increases strongly with T, and that 
this leads in turn to a strong T-dependence of the elastic constants. It emerged 
from this that t he observed anisotropy of sei s mic velocities in the inner core 
fjCreagerl (|l992l ): ISong and HelmbeTgerl (|l993l ): iTromnl (|l993l) ) could be inter- 
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preted in terms of a partial alignment of hep crystallites. Remarkably, the 
calculations required the hep basal plane to be preferentially aligned along 
the Ea rth's rotational axis, which i s the opposite of what had bee n proposed 
earlier ( Stixrude and Cohenl (|l995f ): lsteinle-Neumann et al. f 1999i l). Since the 
correctness of this interpretation depends heavily on the predicted variation of 
c/a with T, it is clearly essential to be confident that the PIC approximation 
can be relied on to give this variation correctly, and one of our aims in this 
paper is to test this point. 



An objective assessment of the errors incurred by the PIC approximation is 
made possible by the fact that ah initio free energies and other thermodynamic 
functions of high-temperature solids can now be calculated with statistical- 
mechanical errors that can be made arbitrarily small. The new methods, de- 
scribed in detail in previous papers (|Alfe et alJ (|200lll2002jl ). are based on the 



DFT calculation of phonon frequencies in the harmonic ap proximation, sup 



pleme nted by the 'thermodynamic integration' technique ([Frenkel and Smit 



(1996)) for calculating anharmonic contributions. For a given material, and 



with a given DFT technique for calculating the electronic total energy as a 
function of ionic positions, we can therefore compute thermodynamic func- 
tions either with the PIC approximation or almost exactly. The differences 
between the two sets of results represent the errors caused by PIC. This is 
what we do in the present work. Since the main difference between PIC and 
the newer, more exact methods is that the latter fully include vibrational cor- 
relations, we shall refer to these in the following as 'vibrationally correlated' 
calculations. 



Two main claims have been made for the PIC method ( Wasserman et al 



f 19961 )): First, that it provides a simple and reasonably accurate way of in- 



cluding the effect of the anharmonicity of lattice vibrations on thermodynamic 
properties; and second, that, even in the absence of anharmonicity, it is com- 
putationally much less demanding than more precise methods. Our assessment 
of these claims for the case of high-p/high-T Fe will suggest that neither is nec- 
essarily true, but that nevertheless the PIC method does yield surprisingly ac- 
curate predictions for many thermodynamic properties. We shall elucidate the 
reasons why PIC is often accurate. As we shall see, however, the present work 
does not support the prediction that the c/a ratio of hep Fe increases strongly 
with temperature, and this casts doubt on the proposed re- interpret at ion of 
the elastic anisotropy of the Earth's inner core. 



The remainder of the paper is organised as follows. In Sec. 2, we describe 
how we have applied PIC to hep Fe and how we have separated the various 
contributions to the free energy and other thermodynamic functions; at the 
end of Sec. 2, we summarise the DFT techniques. In Sec. 3, we describe the 
details of our free-energy calculations and highlight their implications for the 
temperature dependence of the c/a ratio. We then present results for several 
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thermodynamic quantities calculated in the PIC approxima tion, which we 



compa r e with the earlier PIC re s ults o f Cohe n and co-workers flStixrude et al. 



f)l997h : ISteinle-Neumann et all (|200lL 120021): IWasserman et all jl996|)) and 
with the vibrationally correlated results of lAlfe et aL ()2001 ). Discussion and 
conclusions follow in Sec. 4. 



2 Techniques 



2.1 P article-in- cell approximation 



In classical statistical mechanics, the ab initio Helmholtz free energy of a 
vibrating solid containing N ions is: 



*A1 



-k B T In j-^ J dr 1 . . . dr N exp [-/3t/" A i(ri •••r JV )]| 



where U^ii^i • • • r 7v) is the ab initio total electronic free energy of the system 
when the ionic positions are r 1; . . .r^r and A is the thermal wavelength. It 
should be noted that Fai depends on the volume V and temperature T of 
the system (/3 = l/k-sT), and for an hep crystal it also depends on the axial 
c/a ratio, denoted here by q; for the moment we do not indicate explicitly 
the dependence on V, T and q. Two other points should also be recalled. 
First, the quantity C/ai(i"i, • • • , i"Ar) is a free energy, because the electrons are 
treated as being in thermal equilibrium for each set of ionic positions, at a 
temperature equal to the temperature T of the system as a whole. Second, 
the use of classical statistical mechanics, i.e. the neglect of quantu m effects in 



the nu clear motions, is fully justified at Earth's core temperatures (|Alfe et al 
(|200i ). 



The PIC approach consists in replacing Uai in the above formula by the ap- 
proximate form U^{ c , given by: 



N 



U^ c (t u ...r N ) = U peif + ^(Ui 



i=i 



Here, U per f = {7ai(R-i, • • • Rjv) is the ab initio free energy of the system when 
all ions are at their perfect-lattice positions Rj, and 0(uj) is defined to be the 
change of energy {7ai(R<i, . . . Rj + u*, . . . Rat) — t^perf when ion i is displaced 
from its perfect-lattice position to the position R, + u i5 all other ions being 
held fixed at their perfect-lattice positions. With Uai replaced by U^l c in Eqn. 
(1), the PIC approximation for the free energy is: 



F A P i IC = U peii + Nf™ , (3) 
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where the PIC vibrational free energy per atom is: 



/™ 3 = -fc B rin 



A" 3 J du exp(-/30(u)) 



(4) 



The problem is thus reduced to the calculation of a three-dimensional integral, 
which can be performed numerically. 

Our later analysis depends on a clear separation of harmonic and anharmonic 
contributions to f™, which can be accomplished by considering the series 
expansion of 0(u) in powers of ionic displacement u: 

a/3 ' a/3-y ' af3^5 

where the Greek indices a, (3... indicate Cartesian components. The value of 
/^b C when we retain only the quadratic part of 0(u) is the harmonic vibrational 
free energy in the PIC approximation, denoted by f^l? m - This can be expressed 
as: 

/him = k B T £ H^u/k B T) = 3k B T \n(hu/k B T) . (6) 
Here u u are the three PIC vibrational frequencies, given by det \ Mul5 a p — 

(2) 

= 0, with M the ionic mass. The geometric-mean frequency u is defined 
by Eqn. (6). The anharmonic contribution to the free energy is the 

part of not accounted for by f^ m , so that = - f^ m . 

The PIC results reported in Sec. 3 were obtained by computing 0(u) for a 
set of vector displacements u, and fitting the results using the power series 
expansion of Eqn. (5); in practice, for the range of displacements that occur 
with appreciable probability at temperatures below the melting temperature, 
we find that an extremely accurate fit is obtained if we retain only terms up 
to quartic order in u. To show what parameters enter this fit, we write this 
quartic polynomial explicitly for the case of hep symmetry: 



u) ~ \Mul{ul + ul) + M^j + ^ (3) K - 3«H) 



2! 
1 



K^(ul + uy + K^ut(ul + ul) + K^u 



(7) 



Here, u x and u y are perpendicular Cartesian components in the basal plane, 
the x-axis being oriented towards a nearest-neighbour, and u z is the displace- 
ment along the hexagonal axis. The frequencies of harmonic vibration in the 
basal plane and along the hexagonal axis are u) a and u c respectively. This sym- 
metrised polynomial expression may be obtained to a given order by writing 
down all possible terms in a polynomial of that order, and retaining only those 
terms which leave the expression invariant under all the point symmetry oper- 
ations of the chosen cell. In the case of hep symmetry, these include reflection 
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in the hexagonal (z = 0) plane, rotations of ^ about the hexagonal axis, and 
reflection in the y—z plane. With the above form of 0, the geometric-mean 
frequency Q given by 

3 In uj = 2 In u a + In u c . (8) 

The vibrational free energy is obtained essentially exactly by numerical 
evaluation of the integral in Eqn. (4) on a regular grid. In analysing the an- 
harmonic contributions, it is useful to note that can be expanded as a 

power series in temperature: 

/a P nh C a rm = ^ 2 + 0(T 3 ). (9) 

For temperatures below the melting point, we find that only the term in T 2 
is significant, so that in practice the anharmonic contributions are completely 
specified by the coefficient d as a function of volume V and c/a ratio q. The 
following exact expression for the coefficient is readily obtained: 

d = K^/m^lf + K^Jl2(Mu a u c f + K^/SiMuir - (K^f /3(M^) 3 . 

(10) 

The above methods are used to obtain contributions to the free energy at a 
reference value of q, chosen to be close to the T = equilibrium value. We then 
perform further calculations over a range of values of q, to obtain corrections 
to the free energy components, which are parameterised in terms of V and T. 
The details of this parameterisation are given in Sec. 3.1 and 3.2. Equilibrium 
values of q are found by direct minimisation of the total free energy. 



2.2 Ab initio methods 



The DFT electronic-structure techniques used to perform the PIC calcula- 



tions are essentially identical to those described by (jAlfe et al.1 (|200lL 120021 ) ,. 
The exchange- correlation functional is the generalised gradi e nt approxima- 
tion ( GGA) of Perdew and Wang ( Wang and Perdew ( 1991 ): Perdew et al.l 
(|l992n ). The impleme n tation of DFT is the projec tor augmented wave (PAW) 
scheme (|Blochl ( 1994 ): Kresse and Joubert ( 19991 ) ) , with core radii, augmen- 
tation charge radii etc set to the values reported in lAlfe et all (|2000t) . As be- 
fore, all atomic states up to and including 3p-states are treated as core states, 
but the high-pressure response of 3s- and 3p-states, known to be important 
at Earth's-core conditions, is included via an effective pai r potential ; the a c- 
curacy of this procedure has been demonstrated earlier »- 
Thermal excitation of electrons, also important at cor e conditions, is treated 
with t h e usual finite-temperature formulation of DFT (M ermin ()l965l ): lGillar] 
( 19891 ): IWentzcovitch et aD (|l992l ) ). We used a plane-wave cut-off of 300 eV, 
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as in our previous work. AH calc ulations were performed with the VASP code 
f Kresse and Furthmiiller (|l996al lbh). 



The PIC vibrational potential 0(u) (see Eqn. (2)) should in principle be cal- 
culated by displacing a single atom in an infinite crystal. Because PAW cal- 
culations require periodic boundary conditions, we must instead use supercell 
geometry, so that the displaced atom has periodic images. To adhere to the 
PIC scheme, we must therefore ensure convergence of all results with respect 
to the size of the supercell, as described in the following Section. 



3 Results 



3.1 Free energy of perfect lattice 



DFT results for the free energy of the hep Fe perfect latti ce U pPr f(V,q,T) fo r 



the fixed c/a ratio q equal to 1.60 were reported earlier ( Alfe et al. ( 200lf )) 



for V values from 5.2 to 11.4 A 3 /atom at temperatures from 200 to 10 4 K. 
At ea ch T va l ue, th e results were fitted to a Burch-Murnaghan equation of 
state (Poirier (2000)), the parameters of which were then fitted as polynomial 
functions of T. We have repeated these calculations for the present work, and 
as expected the results are virtually identical to those reported earlier. 



To allow for variation of q, we have performed additional calculations of 
C/perf (V, q, T) for q in the range from 1.48 to 1.72, with V going from 5.5 to 10.5 
A 3 /atom, and T going from 2000 to 8000 K. All technical parameters, such 
as the Monkhorst-Pack sampling set, were kept the same as in the q = 1.60 
calculations. To represent the results, we define the deviation AU peT {(V,q,T) 
of the perfect-lattice free energy from its value at a chosen of q, denoted by 
q : AC/perf {V, q, T) = U peT f(V, q, T) - U peT{ (V,q ,T). In the present case, q 
has the fixed value 1.60 throughout. We find that AU pei f(V, q, T) can be very 
accurately represented at all (V, T) by the quadratic form: 

AU peii (q) = a(q-q )(q-q 1 ) • ( n ) 

At each V value, a and q\ can be accurately fitted as linear functions of T, the 
coefficients of which are in turn fitted linearly in V . These fits give a virtually 
perfect representation of the AU peT f(V, q, T) results. Writing a = a + bT, q\ = 
r + sT and with a = + a^V and similarly for b, r and s, the numerical 
values of the parameters are a (0) = 18.39, a (1) = -1.532, 6^ = -3.070 x 10" 4 , 
= 1.64908 x 10 5 , r<°> = 1.63, = -7.909 x 10~ 3 , = -7.902 x 10~ 6 
and s^ 1 ) = 1.817 x 10~ 6 . Units of coefficients are such that a is in eV, T in K 
and V in A 3 . 
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The equilibrium value g of the perfect lattice, denoted by g eq , is obtained by 
minimising AU peT { (V, q, T) with respect to q. At T = 0, we find a very weak 
dependence of q cq on V, going from 1.590 at 7.0 A 3 /atom to 1.578 at 10.0 
A 3 /atom. The T-dependence of AU peT f(q) is also very weak, with a decreasing 
linearly from 7.4 to 6.2 eV, and gi increasing from 1.59 and 1.62 as temperature 
varies from 2000 to 8000 K. If we now use the results to predict g eq for the 
(hypothetical) high-T perfect lattice, we find a variation of at most 0.04 at a 
volume of 10.0 A 3 . The insignificant V- and T-dependence of g eq is noteworthy, 
because it implies that any significant variation of q eq at high temperatures 
can come only from lattice vibrations, to which we now turn. 



3.2 Vibrational free energy 



The basal-plane and axial frequencies u a and uj c , and the four anharmonic 
vibrational coefficients K^L an d Kf^ were calculated as follows. 

For given values of q and V, a supercell was constructed in which one of the 
atoms, the 'walker' is displaced from its equilibrium position. The 'walker' is 
given a series of equally spaced displacements in the three directions u y , u z 
and ^(u x + u z ), the maximum displacements r max in each direction being 
chosen so that the Boltzmann factor exp [— /90(r max )] « 0.1 for the maximum 
temperatures of interest. Directions were selected to allow the terms in 0(r) 
(Eqn. (7)) to be determined as simply as possible by least squares fitting. 
Size convergence in terms of supercells is essential to ensure that the walker 
cannot 'see' its image. We discuss the results of a simple convergence test 
below. All the calculations described here are performed at a fixed electronic 
temperature of 6000 K, since statistical mechanical considerations will dom- 
inate the temperature dependence of /^b C - At this temperature, vibrational 
coefficients were fully converged with respect to fc- point sampling at each cell 
size. Monkhorst-Pack ( Monkhorst and Pack (197^)) sampling was used, with 



9x9x5 /c-points for the 8 atom cell, 7x7x5 /c-points for 16 atoms and 
3x3x3 for 36. Calculations were performed over a range of volumes from 
5.5 A 3 /atom to 11.5 A 3 /atom with a fixed axial ratio of go — 1-60 as before. 
Using 5th order polynomials to obtain a high quality fit, the quantities if", 
K^ x , and Inu were then parameterised in terms of volume. Finally, 
a correction due to relaxation of q is added to the harmonic free energy, as for 
that of the perfect lattice, by performing a parameterisation of Q in terms of 
q, which we describe below. Parameterisation of higher order K coefficients in 
terms of q yielded only negligible corrections to the p-V curve. 

As noted in Sec. 2.1, the harmonic free energy is completely determined by 
the geometric-mean frequency u (Eqns. (6) and (8)), and we consider first our 
PIC results for uo as a function of V for the fixed c/a ratio go — 1-6. We report 
in Fig. 1 our ui results from calculations using supercells of 8, 16 and 36 atoms, 
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compared with the Co result s from our earlier calculations of the full phonon 
spectrum ( Alfe et al. ( 200lh ) . We note two important points. First, the PIC Co 
results are almost independent of supercell size, so that the calculations appear 
to be fully converged with respect to size effects. Second, the PIC Co differs 
from the full-phonon Co only by a volume-independent shift of In Co. Since such a 
^/-independent shift cannot affect most thermodynamic quantities, we expect 
the PIC calculations to agree rather closely with full-phonon calculations. 



We now consider the dependence of Co on c/ a ratio q, defining the logarithmic 
deviation Alu>(V, q) from its value at go by lnu)(V, q) = \nCo(V, q ) + Ai^Co{V, q). 
In terms of this, the harmonic free energy can be written as: 



fZ^JV, q, T) = 3k B T \n(hio/k B T) = f^jV, q Q , T) + 3A; B TA L ^ . (12) 

Calculations of Co(q) were made for q- values going from 1.48 to 1.72 at a se- 
ries of volumes. We find that Co decreases with increasing q. Not surprisingly, 
increase of q at constant V yields an increase of the basal-plane frequency uo a 
and a decrease of axial frequency uo c . The decrease of uo c succeeds in outweigh- 
ing the increase of uo a in the formula for Co (Eqn. (8)), but nevertheless the 
resulting decrease of Co is a fairly marginal effect. Over the g-range studied, we 
found that the quadratic form A L a>(g) = (3(q — qo)(q — q 2 ) gives an accurate fit, 
and this fit was performed at volumes of 5.5, 7.0 and 10.5 A 3 /atom. The In- 
dependence of the resulting (3 and q 2 coefficients was obtained by a quadratic 
fit. Writing (3 = J2i=o P^V 1 and similarly for q 2 , the values of the parameters 
are = -2.083, /?W = 0.375, /3® = -0.0264, q { 2 0) = 1.990, q£ ] = -0.141 
and q 2 = 8.79 x 10~ 3 . Units of parameters are such that Co is in rad s _1 and 
V is in A 3 . 



From the values of the four anharmonic coefficients at each V and q, we 
calculate also the anharmonic vibrational free energy /^L m . As explained 
in Sec. 2.1, the most accurate way of doing this is by explicit calculation of 
/^ b c (Eqn. (4)) by numerical evaluation of the integral over displacement u, 
from which /^ha rm is obtained as the difference — /^^. Alternatively, 
we can use Eqn. (10) to calculate the d coefficient of the leading term in 
the temperature expansion of /^hi rm (see Eqn. (9)). In practice, we find that 
the two methods yield almost indistinguishable results, and everything that 
follows is based on evaluation of the d coefficient. 



We report in Fig. 2 the volume dependence of d for q = 1.60, compared with 
the corresponding anharmonicity co efficient obtained from the vibrationally 
correlated results of lAlfe et al.l ([20011 ). The coefficients of this fit are as follows: 
= -6.556 x 10" 8 , d^ = 3.399 x 10" 8 , = -6.487 x 10~ 9 , d® = 
5.400 x 10~ 10 , d^ = -1.639 x 10" 11 , where d(v) = Y.L n d®V\ d be ing in 
eV K -2 and V in A 3 . Corresponding data is given in lAlfe et al.l ([20011 ). The 
striking and important feature of this comparison is that the PIC anharmonic 
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free energy is positive, whereas the results of lAlfe et al.l (|200l[ ) show that in 
reality it is negative. This means that for high-p hep Fe, the PIC approximation 
gives a completely incorrect account of anharmonicity, the electronic structure 
techniques in both works being identical. However, we note that f^harm * s * n 
any case very small. For example, at V = 7 A 3 /atom, T = 6000 K, /^harra ~ 
15 meV/atom. Because it is so small, we ignore its dependence on q in the 
following, and use a polynomial fit to the d(V) results for q = 1.60. 



3.3 Temperature dependence ofc/a ratio 



The equilibrium value q cq of the c/a ratio q at given V and T is obtained by 
minimising the total free energy Fj^ c with respect to q. The variation of q with 
T comes from the T- and g-dependence of the perfect-lattice free energy , 
for which we presented results in Sec. 3.1, and from the g-dependence of the 
mean vibrational frequency u (Eqn. (6)). (Since we neglect the g-dependence of 
the anharmonic component of free energy, anharmonicity does not contribute 
to the T-dependence of g eq here.) 

Our results for g eq at three different volumes are re ported in Fig. 3, where 



we compare them with the earlier PIC predictions of ISteinle-Neumann et al. 



(200ll 12002). The present results are very different from the earlier ones. At 



all volumes, we find only a very weak increase of q w ith T, which is between 



5 and 10 times smaller than the variation reported in ISteinle- Neumann et al 



( 20011 ). We note that this gross discrepancy can come only from a difference 
in the g-dependence of the harmonic mean frequency Q. The reason for this is 
that the roughly linear dependence of g cq on T seen in both our results and 



those of ISteinle- Neumann et al.l ((2001) cannot originate either from thermal 



electronic excitations or from anharmonicity, since both of these free-energy 
components vary as T 2 . The cause of the discrepancy, and its implications for 
understanding the elastic anisotropy of the inner core, are discussed further 
in Sec. 4. 



3.4 Thermodynamic functions 



All thermodynamic functions can be calculated by taking appropriate deriva- 
tives of the PIC free energy , with its perfect-lattice, harmonic and an- 
harmonic components parameterised as described above. All the results to 
be presented include the dependence of the equilibrium c/a ratio g eq on V 
an d T. Our PIC pr e dictio ns will be compared with the earlier PIC results 
of Wasserman et al. ( 1996| ) and with the vibrationally correlated results of 



Alfe et al.l (|200ll ). 
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We begin by considering the thermal pressure Ap, which is one of the most 
basic quantities in interpreting the properties of the inner core. This is defined 
as the difference between the pressure in the material at a given V and T and 
the pressure at the same volume but at zero temperature: 



Ap = -(dF^/dV) T + (dF^/dV) T=0 . (13) 
Fig. 4 shows a com parison of our PIC results for Ap with the vibrationally 



correlated results of lAlfe et al.l (|2001[ ) on isotherms at T = 2000, 4000 and 



6000 K. The rather close agreement shows that the PIC approximation gives an 
accurate account of the thermal pressure. This is expected from our results for 
u, since in the harmonic approximation Ap is given by the electronic thermal 
pressure, which is exactly the same in both calculations, plus a vibrational 
contribution equal to —3k^Td In ui/dV. We have seen that for hep Fe, the 
PIC value of In a) differs from the exact value by an almost constant offset, 
which has no effect on the derivative of In u with respect to V. Small differences 
between the PIC and vibrationally correlated thermal pressure may be due to 
small variations in the offset between the values of In a). 

The situation is similar for the thermal energy AE, defined as the difference 
between the internal energy E at a given V and T and the internal energy at 
the same V but zero temperature: AE = E(V,T) — E(V,0). The electronic 
part of the thermal energy is identical in both PIC and vibrationally corre- 
lated calculations, and the harmonic vibrational energy is exactly 3ksT / atom 
in both cases, so that any difference arises only from the small anharmonic 
contribution. Results for the thermal energy are not presented here. 

The good accuracy of PIC for thermal pressure and thermal energy explains 
why it has also been fou nd to give a good account of shock measurements 



( Wasserman et al.l (J1996J)). A conventional sequence of shock experiments on 



samples in the same initial state generates a path through thermodynamic 
state-space known as the Hugoniot. On this path, the pressure p, vo lume V 



and in ternal energy E are related by the Rankine- Hugoniot formula ([Poirier 
(|200nh ): 

~(p + Po)(V -V) = E-Eo, (14) 

where p , V and E refer to the initial state, which usually corresponds to 
ambient conditions. For given V, the Hugoniot p and temperature T can be de- 
termined by going to the calculated p(V, T) and E(V, T) functions and seeking 
the value of T for which p and E satisfy eqn (14). Fig. 5 shows our calculated 
PIC p(V) on the Hugoniot co mpared wi t h the experimental results and with 
the fully correlated results of lAlfe et al. (|200lh . The almost exact agreement 



between the three sets of results confirms the excellence of PIC for this par- 
ticular purpose. Similar comparisons are shown for the T(p) relation on the 
Hugoniot in Fig. 6. Here there is a discrepancy, which corresponds to a dif- 
ference in internal energy between the two theoretical models. This difference 
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is fully accounted for by the neglect of temperature dependence of u> in the 
present work. 



We conclude the presentation of results by examining three thermodynamic 
quantities for which our comparis ons between vibrat i onally cor related calcula- 
tions a nd the earlier PIC results of Washerman et all (Il99fih and Stixrude et al 



(1997|) revealed significant differences ([Alfe et al.l ([20011 )) 



Fig. 7 compares the present PIC results for thermal expansivity a with earlier 
calculations. At 2000 K, the present results are in very close agreement with 
the vibrationally correlated results, and in strong disagreement with the earlier 
PIC results at low pressures. At higher T, the present results fall somewhat 
below both previous sets of results. The product aK T (Fig 8) of expansivity 
and isothermal bulk modulus is important in high-pressure work, because it 
can sometimes be assumed to be independent of p and T over a wide range 
of conditions. The vibrationally correlated calculations showed that for high- 
p/high-T hep Fe constancy with p is a good approximation, but constancy 
with T is not. The present PIC results show a reduction of aKnp with respect 
to vibrationally correlated and earlier PIC results of, at most, around 15% at 
6000K. 

The thermodynamic Griineisen parameter 7 = V(dp/dE)v is particularly im- 
portant, because it relates thermal pressure to thermal energy, and assump- 
tions about 7 are often used in reducing shock data from Hugoniot to isotherm. 
The large differences between the earlier PIC results an d the vibra t ionall y- 
correlated results for 7 are therefore a cause for concern (|Alfe et al.1 jiioi ). 



The present PIC results (Fig. 9) agree quite closely with the vibrationally- 
correlated results, and this suggests that the cause of the earlier disagreement 
was not the PIC approximation itself. 



4 Discussion 



Our comparisons with more exact calculations have shown that the PIC ap- 
proximation gives very good results for a range of important thermodynamic 
properties of hep Fe at Earth's core conditions. A noteworthy example of this 
is that piV) on the Hugoniot agrees almost perfectly with the more exact 
results. Our analysis of the free energy into different components makes clear 
why PIC is so good. The perfect-lattice component is exactly the same in the 
two approaches. For the harmonic component, the sole requirement for good 
results is that the logarithmic derivative d In uo/d In V of the geometric mean 
frequency uj be correct. But we have seen that for hep Fe the PIC uj differs 
from the uj given by calculation of the full phonon spectrum by an almost 
constant factor over a wide range of volumes, so that this requirement is sat- 
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isfied. The basic reason for this is that the phonon dispersion relati ons of hep 
Fe sca le by a uniform factor with changing volume (see Fig. 3 of lAlfe et al. 



(|200l ). " ? inally, the anharmonic component of free energy is very small, and 
has only a very minor effect on most thermodynamic functions. The reliabil- 
ity of PIC actually requires that anharmonic effects be small, since we have 
shown that the PIC approximation misrepresents these effects in predicting 
the wrong sign of the anharmonic free energy. 

Surprisingly, even though PIC appears to be so good, we fin d im portant dis- 

crepan cies with the earlier PIC results of lStixrude et al.l ( 1997 ) and lWasserman et al 
(1996). In particular, our calculations of the thermod ynamic Gruneisen param- 
eter agree much more closely with the calculations of lAlfe et al.l ((20011 1 . These 
discrepancies are clearly not due to PIC itself, b ut must come from othe r 



technical differences. We note that in the work of Wasserman et al. (1996) 



the PIC calculations actually employed a tight-binding representation of the 
total energy function, the parameters in the tight-binding model being fitted 
to ab initio calculations. Conceivably, the tight-binding fit might have led to 
errors. 



Even more surprising is that the strong increase with temperature of the axial 
c/a ratio predicted by recent PI C calculations is not reproduc ed at all by the 
present PIC work. According to ISteinle- Neumann et al.l ( 2001 ). at the atomic 
volume of 7.11 A 3 , c/a increases from 1.63 to 1.75 as T goes from 2000 to 
8000 K, whereas in the present PIC calculations at the similar volume of 
7.0 A 3 , c/a increases only from 1.594 to 1.610 over the same temperature 
range. The correctness o f the weak cj a variation found her e is suppor t ed by 
prelim inary calculations ([Alfel ( 2002| )) using the techniques of Alfe et al. (2001, 
20021) . The reasons for this discrepancy are completely unclear at present. The 



discrepancy has major implications for our understanding of the Earth's inner 
core, because the recently proposed re-interpretation of the elastic anisotropy 
of the inner core appears to depend crucially on a strong T variation of c/a. 
We believe it is highly desirable that this question be re-examined by other 
research groups. 

Our work sheds light on the usefulness of the PIC approximation. Since we 
have seen that PIC cannot be relied on for anharmonic contributions, it should 
be regarded as a way of calculating the geometric-mean harmonic frequency 
u. But PIC requires ab initio calculations on a periodic system in which a 
single atom is displaced from its perfect-lattice site. These are precisely the 
same calculations that are performed in order to obtain the force-consta nt 



20021) ). 



matrix used to compute the full phonon spectrum (|Alfe et al.l (|2001 , 
For an ab initio method that yields forces on all ions - which includes the 
pseudopotential and PAW implementations of DFT, among others - the net 
result of the PIC procedure is to discard all the information contained in the 
ionic forces, retaining only the variation of total energy with displacement. 
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This suggests that it may be better to use the force information to compute 
the force-constant matrix and hence the full phonon spectrum, rather than 
adopting the PIC approach. 
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Fig. 1. The variation with atomic volume V of In a) (a) is the geometric mean vibra- 
tional frequency) from PIC calculation on periodic cells containing 8 (so lid curve), 
16 (da shed) and 36 (dotted) atoms. The vibrationally correlated results of lAlfe et al 
(2001) are given by the lighter curve. 
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Fig. 3. Equilibrium axial ratio q, for lSteinle-Neumann et all ( 200lh (heavy curves) at 
atomic volumes of 6.8lA 3 (solid curve), 7.1lA 3 (dashed curve) and 7.4lA 3 (dotted 
curve), and for the current work (light curves) at 5.5A 3 (solid curve) 7.5A 3 (dashed 
curve) and 10. OA 3 (dotted curve) 
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Fi g. 4. The t o tal th ermal pressure on isotherms in this work (light curves) and that 
of klfe et al.1 (fcOQlh (heavy curves) at 2000 K (solid), 4000 K (dashed) and 6000 K 



(dotted) function of atomic volume. 
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Fig. 6. The temp erature-pres s ure H ugoniot. Heavy and light curves correspond 

2001) respectively; black circles show the exper- 
(1993) and empty circles are estimates due to 



to this work and Alfe et all 
imental results of Yoo et alJ 
Brown and McQueen! ()l98fil ). 
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Fig. 7. Thermal expansivity on isotherms at 2000 K (solid) an d 6000 K ( dotted 
curves). Heavy, medium and li ght curves corre s pond to th is work, Alfe et alj 
and the earlier PIC results of 
respectively. 



j;nt curves correspond to this work, Alte et aJ 
IStixrude et al. (|l997l ) and IWasserman et al 
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Fig. 8. The product aKx on isotherms at 2000 K (solid) and 6000 K 



Heavy, medium and li ght curves corre s pond to thi s work, Alfe et al 



earlier PIC results of iStixrude et al 
tively. 



and IWasserman et al 



dotted curves) . 
and the 
respec- 
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Fig. 9. The Griineisen parameter on isotherms at 2000 K (solid) a nd 6000 K (dotted 



curves). Heavy, medium and li ght curves corre s pond to th is work, Alfe et al.j 
and the earlier PIC results of 
respectively. 



yit curves correspond to this work, Alte et aJ 
IStixrude et al. (|l997l ) and IWasserman et al 
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